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Abstract The penalization method is used to take account of obstacles in 
a tokamak, such as the limiter. We study a non linear hyperbolic system 
modelling the plasma transport in the area close to the wall. A penalization 
which cuts the transport term of the momentum is studied. We show numer- 
ically that this penalization creates a Dirac measure at the plasma-limiter 
interface which prevents us from defining the transport term in the usual 
sense. Hence, a new penalty method is proposed for this hyperbolic system 
and numerical tests reveal an optimal convergence rate without any spurious 
boundary layer. 
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1 Introduction 



A tokamak is a machine to study plasmas and the fusion reaction. The plasma 
at high temperature (10 8 K) is confined in a toroidal chamber thanks to a 
magnetic field. One of the main goals is to perform controlled fusion with 
enough efficiency to be a reliable source of energy. But, since the magnetic 
confinement is not perfect, the plasma is in contact with the wall. In order 
to preserve the integrity of the wall and to limit the pollution of the plasma, 
it is crucial to control these interactions. 

We study, using a fluid approximation of the plasma, a simplified system 
of equations governing the plasma transport in the scrape-off layer, parallel 
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to the magnetic field lines. In this paper, after a numerical study of the pena- 
lization introduced by Isoardi et al. [9] , we modify the boundary conditions to 
ensure the well-posedness of the hyperbolic system and we propose another 
penalty method which seems to be free of boundary layer. 



2 The model hyperbolic problem 

In this paper, we consider a very simple model taking only into account the 
transport in the direction parallel to the magnetic field lines, (see for example 
[9], [13]). It is a one dimensional 2x2 hyperbolic system of conservation laws 
for the particle density N and the particle flux r, which reads : 



Here, the boundaries of the domain x — L and x = —L correspond to the 
"limiters", which are material obstacles for the fluid (see Fig. 1). In the right- 
hand side, S is a source term. 

There is a difficulty with the choice of the boundary conditions for the 
system ([!]). From physical arguments, it follows that the domain (namely the 
scrape-off layer) is basically divided into two regions [13] : 

• One region far from the limiter, the pre-sheath, where the plasma is neutral 
and the Mach number M = r/N of the plasma satisfies \M\ < 1. 

• One region next to the limiter (in a thin layer called the sheath area, whose 
typical thickness is of the order of 10 _5 m), where the electroneutrality 
hypothesis does not hold and we have \M\ > 1. More precisely M > 1 
close to x = L and M < —1 close to the boundary x = —L. 

It could seem natural to prescribe M — 1 (resp. M = —1) as a boundary 
condition at x = L (resp. x = —L) for the system, since the physical argu- 
ments imply that M = ±1 very close to the obstacle (Bohm criterion). These 
are exactly the boundary conditions which are chosen in [9]. However, in that 
case, as the eigenvalues of the Jacobian of the flux function are M — 1 and 
M + 1, it follows that at the plasma limiter interface one eigenvalue is (the 
boundary is characteristic) and the other one is outgoing (it is also true at 
x = —L), and clearly the problem does not satisfy the usual sufficient condi- 
tions for well posedness, see [8] , [IT] , [3] : the number of boundary conditions 
(= 1) is not equal to the number of incoming eigenvalues (= 0). 

In order to test our penalty approach with a well-defined hyperbolic bound- 
ary value problem, in section [3] we slightly modify the boundary conditions 
of the paper [3], and impose M=l — eonx = L and M = — 1 + e on x = —L 




Initial conditions : N(Q, .) = A?o and T(0, .) = r Q 
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with a fixed e > 0, which leads to a well-posed hyperbolic problem. In our 
numerical simulations we use e = 0.1. 




Fig. 1 Schematic representation of the scrape-off layer. The x-axis corresponds to the 
curvilinear coordinate along a magnetic line close to the wall of the tokamak. 



The numerical tests presented below, use a finite volume scheme with 
a second order extension : MUSCL reconstruction with the minmod slope 
limiter and the Heun scheme which is a second order Runge-Kutta TVD 
time discretization. The finite volume scheme is the VFRoe using the non 
conservative variables for the linearized Riemann solver [7j; here, the non 
conservatives variables are N and M. To avoid stability issues, the penalized 
terms are treated implicitly for the time discretization. 



3 Study of penalty methods 
3.1 A first penalty method 

The following penalty approach has been proposed by Isoardi et al. [3]. Let's 
X be the characteristic function of the limiter, i.e. x( x ) — 1 if x is in the 
limiter, and x( x ) = elsewhere, and r\ the penalization parameter. The 
penalized system is given by : 

{d t N + d x r+^N=(l-x)S N inR+xK 
r) 
d t r + (1 - x )d x (j^+ N )+^ r - M oN) = (1 - x)Sr (2) 
Initial conditions : N(0, .) = 7V and f(0, .) = T 

Mo is a function such that, at the plasma-limiter interface we have \Mq\ = 1. 
Here, the two components of the unknown are penalized although there is 
no incoming wave. At least formally, N is forced to converge to inside the 
limiter when r\ tends to 0. 
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The flux of the second equation is cut inside of the limiter, and this causes 
some troubles from the mathematical point of view. Indeed, this is an hyper- 
bolic system with discontinuous coefficients and the meaning of the term 

(1-X)d x (jf+ N ) 

is not clear because it can involve the product of a measure with a discon- 
tinuous function which has no distributional sense. As a consequence and 
as a confirmation of this fact, our numerical tests show the existence of a 
strong singularity at the interface for the numerical discrete solution. Con- 
cerning the interpretation of this numerical singularity, it could happen (but 
we don't have any rigorous proof and this is just an open question) that this 
system admits generalized solutions in the spirit of Bouchut-James [4] (see 
also Poupaud-Rascle [TO], or Fornet-Gues [3]) such as measure- valued solu- 
tions, which can for example exhibit a Dirac measure at the interface, and 
this generalized solution could be selected by the numerical approximation 
process. 

For the numerical test, we choose Sjsr and Sr so that the following functions 
define a solution of the boundary value problem : 

N ^ X) = exp (oie^i)) rM = sin © exp (oiiri)) 

These test solutions are regular (at least inside the plasma area) and has no 
singularity at the plasma-limiter interface. In the Fig. [2j we observe that a 




Fig. 2 M versus x with rj = 10 — 3 , a mesh of J = 1280 cells using the penalization of 
Isoardi et al. [9]. The computations are stopped when maxigj-^ j}(|Af™|) > 10, which 
corresponds to the time : t = 0.008822. The computational domain was [0, 0.5] and L = 0.4 
(plasma-limiter interface). At x = 0, we impose a symmetry condition. 



3 Study of penalty methods 



5 



peak appears very quickly, then \M"\ become very large (about 10 8 ) in a 
few points. The same computations are made for two more refined meshes 
(respectively for 2560 and 10240 cells) and we observe that the peak is nearer 
and nearer to the plasma limiter interface, when the resolution increases. 
Besides, when the mesh step decreases, the peak appears earlier and earlier. 
We stop the computations when maXj e {i,...n(|M™|) > 10 but similar results 
are obtained when the stop criterion is max^g^ > 100. This leads 
one to believe that, if the solution converges to a generalized solution of the 
continuous problem, then this generalized solution must have a singularity 
supported by the interface (that could be a Dirac measure for example). We 
notice that the presence of a Dirac measure at the interface is not only a 
theoretical issue since it has been observed numerically and that the Dirac 
measure destabilizes numerical schemes. In the following section, we propose 
a modification of the boundary value problem to obtain a well-posed version. 



3.2 A new penalty method for the modified boundary 
conditions 

After the modifications proposed in section [2] the well-posed initial boundary 
value problem reads : 



For this problem, the boundary is not characteristic, and the boundary condi- 
tions are maximally dissipative. Hence, for compatible initial data, the prob- 
lem has a unique local in time solution, which is regular enough : at least C 1 
is sufficient to perform the asymptotic analysis; see e.g. [3J, [T^] - 

To penalize ([3]), we use a method developed in the semi-linear case by 
Fornet and Gues [BJ. In order to have an homogeneous Dirichlet boundary 
condition for the theoretical study, the system is reformulated with the un- 
knowns u = ln(iV) and v = r/N — Mq. Although our system is quasi-linear 
(and not semi-linear), the method can be extended to this case. An interesting 
feature of the method is that it yields to a convergence result without gener- 
ation of a boundary layer inside the limiter. Up to now, we don't know if this 
method can be extended to more general quasi-linear first order hyperbolic 
system with maximally dissipative conditions. 

We assume that Mq is a constant such that < Mq < 1. We denote by \ 
the characteristic function associated to the limiter, i.e. xi x ) = 1 if the point 
x is in the limiter. 




M (., —L) = -1 + e and M(.,L) = 1 - e 
k JV(0,.) =N and T(0, .) = T 
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The new penalized problem reads : 

(dtN + a x r = s N 

U r+9 *(^ +iV ) + f {w - N ) =Sr mRtxR (4) 

[ JV(0, .) = JV and r(0, .) = r 

The formal asymptotic expansion of a continuous solution to Q with the 
BKW (Brillouin-Kramers-Wentzel) method does not contain any boundary 
layer term [1] and this suggests strongly that there is no boundary layer at 
all in the solution. Notice that the penalization is incomplete: only one field 
is penalized, which is natural since there is only one boundary condition. 

For the numerical tests, we use a regular solution: 

N(f > x) = exp ( o.ieft + i) ) nt ' x) = MoSin O exp ( o-ieft + i) ) 

and Siy,Sr are well chosen. The spatial domain is [0,0.5] with a symmetry 
condition at x = and the limiter set corresponds to x <E [0.4, 0.5]. 



Thick black : N(1 .x). Gray : Gamma(1 ,*), Narrow black : M(1 ,x), ela=1e-1 




Fig. 3 Plot of N, r and M as functions of x (at t = 1) with the penalty method free of 
boundary layer for 77 = 0.1. The continuous lines represent the numerical solutions whereas 
the dashed lines corresponds to the exact solution of the hyperbolic limit problem (tj — > 0). 
The limiter corresponds to the area x S [0.4,0.5]. For smaller values of r/, for instance for 
r] = 10~ 5 , the plot is almost the same as the plot of the exact solution (dotted lines). 



We analyze the convergence when the penalization parameter r\ tends to 
using a uniform spatial mesh of step Sx = 10~ 5 . We calculate the error 
in L 1 norm for N, d x N, r and d x T. The goal is to confirm numerically the 
absence of boundary layer with an optimal rate of convergence as 0(rj). 
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One of the main difficulties for the implementation of the penalization, 
is the choice of a boundary condition at x = 0.5 which is necessary for the 
numerical scheme. As only r is penalized, we need a transparent boundary 
condition for N. For the numerical tests, the boundary condition comes from 
the asymptotic expansion up to the first order of the BKW analysis. We 
carry out the computations up to t = 1 with an adaptive time step so that 
the CFL condition is always satisfied. The results are plotted in Fig. [3j In 
Fig. El we observe that the optimal rate of convergence 0(rj) is reached for 
the Lr norms, even for the derivatives. This gives a numerical evidence of the 
absence of boundary layer. The same numerical results in 0(r]) are obtained 
if the penalty term in Q is replaced by ^ — M ), see [2]- 
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When the parameter e = 0.01, i.e. close to a characteristic boundary, the 
computations show that, for 77 sufficiently small, 77 < 0(e), the convergence 
results are similiar; see details in [T]. 
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